{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Iterative Template Generation\n",
    "\n",
    "This notebook exemplifies one way in which a template mesh atlas can be generated from a collection of segmented binary images. Each binary image of a mouse femur is downsampled to reduce pixel density prior to applying marching cubes to generate a mesh from the binary image. One arbitrary mesh is selected as the template and then registered to and resampled from each original mesh to get a full set of meshes with correspondence points. The meshes are then groupwise registered via procrustes alignment and the mean mesh is taken as the new template. This process is repeated for a fixed number of iterations to get a template mesh atlas that represents the average case of all meshes.\n",
    "\n",
    "This pipeline uses the [ITKShape](https://github.com/slicersalt/ITKShape) module for shape analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "!{sys.executable} -m pip install numpy itk itk-shape itkwidgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import time\n",
    "\n",
    "import itk\n",
    "import numpy as np\n",
    "from itkwidgets import view, checkerboard\n",
    "\n",
    "module_path = os.path.abspath(os.path.join('..'))\n",
    "\n",
    "if module_path not in sys.path:\n",
    "    sys.path.append(module_path)\n",
    "\n",
    "from src.hasi.hasi import align\n",
    "from src.hasi.hasi.pointsetentropyregistrar import PointSetEntropyRegistrar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read images\n",
    "\n",
    "Input images represent the results of automatic binary segmentations of mouse femur data. Each image contains only the femur object and potentially represents a different region in space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "IMAGE_FOLDER = 'Input/femurs/'\n",
    "DENSE_MESH_OUTPUT_FOLDER = 'Output/femurs/'\n",
    "TEMPLATE_OUTPUT_FOLDER = 'Output/templates/'\n",
    "MEAN_OUTPUT_FOLDER = 'Output/mean/'\n",
    "\n",
    "for folder in [IMAGE_FOLDER, DENSE_MESH_OUTPUT_FOLDER, TEMPLATE_OUTPUT_FOLDER, MEAN_OUTPUT_FOLDER]:\n",
    "    os.makedirs(folder, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Input/femurs\\\\901-R-femur-label.nrrd', 'Input/femurs\\\\902-R-femur-label.nrrd', 'Input/femurs\\\\906-R-femur-label.nrrd', 'Input/femurs\\\\907-R-femur-label.nrrd', 'Input/femurs\\\\908-R-femur-label.nrrd', 'Input/femurs\\\\915-R-femur-label.nrrd', 'Input/femurs\\\\916-R-femur-label.nrrd', 'Input/femurs\\\\917-R-femur-label.nrrd', 'Input/femurs\\\\918-R-femur-label.nrrd', 'Input/femurs\\\\F9-3wk-01-R-femur-label.nrrd', 'Input/femurs\\\\F9-3wk-02-R-femur-label.nrrd', 'Input/femurs\\\\F9-3wk-03-R-femur-label.nrrd', 'Input/femurs\\\\F9-8wk-01-R-femur-label.nrrd', 'Input/femurs\\\\F9-8wk-02-R-femur-label.nrrd']\n"
     ]
    }
   ],
   "source": [
    "# Get healthy femur segmentation binary images at \n",
    "# https://data.kitware.com/#collection/5dcc6691e3566bda4b802172/folder/5e0b8d6baf2e2eed35c326f7\n",
    "\n",
    "input_paths = glob.glob(IMAGE_FOLDER + '*-R-*')\n",
    "assert(len(input_paths) == 14)\n",
    "\n",
    "print(input_paths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "MESH_FILENAMES = [os.path.basename(file).replace('.nrrd','.obj') for file in input_paths]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "images = list()\n",
    "for path in input_paths:\n",
    "    images.append(itk.imread(path, itk.UC))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Paste images into same space\n",
    "\n",
    "Here we standardize physical space across the femur images. This is primarily intended to assist in viewing convenience with `itkwidgets` which expects a standard viewing region, but could also be helpful to standardize output from subsequent image downsampling and mesh conversion operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "images = align.paste_to_common_space(images)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d6416cb02cf2417a89138e9344b3dfbd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[], gradient_opacity=0.22, point_sets=[], rendered_image=<itk.itkImagePython.itkImageUC3; pr…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view(images[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Downsample images\n",
    "\n",
    "The marching cubes algorithm returns a mesh with vertex density related to the pixel density of the original image. Here we downsample each image twice, once to get a \"dense\" image retaining most information density and a second time to get a \"sparse\" image more easily applied to get correspondence points.\n",
    "\n",
    "Meshes generated later from the dense images have approximately 600,000 vertices each while meshes generated from the sparse images have approximately 4,000 vertices. We will use the dense meshes to sample feature information and iteratively refine a the atlas to generalize the shape population. We can select a single sparse mesh to act as the initial atlas or carry out iterative refinement on multiple sparse meshes and compare to determine which result \"best\" reflects the population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "SPARSE_DOWNSAMPLE_RATIO = 14\n",
    "DENSE_DOWNSAMPLE_RATIO = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "itkSize3 ([639, 477, 519])\n",
      "itkSize3 ([91, 68, 74])\n"
     ]
    }
   ],
   "source": [
    "sparse_downsampled_images = align.downsample_images(images,SPARSE_DOWNSAMPLE_RATIO)\n",
    "dense_downsampled_images = align.downsample_images(images,DENSE_DOWNSAMPLE_RATIO)\n",
    "\n",
    "print(itk.size(dense_downsampled_images[0]))\n",
    "print(itk.size(sparse_downsampled_images[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6f25616a386b4871b8e0a328105932ea",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[], gradient_opacity=0.22, point_sets=[], rendered_image=<itk.itkImagePython.itkImageUC3; pr…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view(dense_downsampled_images[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Meshes\n",
    "The `itk.BinaryMask3DMeshSource` class makes use of the Marching Cubes algorithm to generate a mesh from a given object. Each binary image here uses the value '1' to indicate the femur is present at a pixel and '0' to indicate the femur is not present. Marching Cubes rapidly fills the femur space and generates surfaces at pixel region boundaries.\n",
    "\n",
    "Note that it may be useful to visually examine intermediate results. In the case where a mesh is not well aligned with others it is useful to correct the transformation in an external mesh editor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Here each pixel interior to the femur has value \"1\" and exterior has value \"0\".\n",
    "# This may change for a different segmentation image.\n",
    "FEMUR_OBJECT_PIXEL_VALUE = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average dense mesh points: 215237.7857142857\n"
     ]
    }
   ],
   "source": [
    "dense_meshes = align.binary_image_list_to_meshes(dense_downsampled_images,\n",
    "                                           object_pixel_value=FEMUR_OBJECT_PIXEL_VALUE)\n",
    "\n",
    "# Expect ~200K vertices\n",
    "print('Average dense mesh points: ' +\n",
    "      str(sum([mesh.GetNumberOfPoints() for mesh in dense_meshes]) / len(dense_meshes)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average sparse mesh points: 4371.5\n"
     ]
    }
   ],
   "source": [
    "sparse_meshes = align.binary_image_list_to_meshes(sparse_downsampled_images,\n",
    "                                            object_pixel_value=FEMUR_OBJECT_PIXEL_VALUE)\n",
    "\n",
    "# Expect ~200K vertices\n",
    "print('Average sparse mesh points: ' +\n",
    "      str(sum([mesh.GetNumberOfPoints() for mesh in sparse_meshes]) / len(sparse_meshes)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Write out each mesh to disk\n",
    "for i in range(0,len(dense_meshes)):\n",
    "    itk.meshwrite(dense_meshes[i], f'{DENSE_MESH_OUTPUT_FOLDER}{MESH_FILENAMES[i]}')\n",
    "\n",
    "for i in range(0,len(sparse_meshes)):\n",
    "    itk.meshwrite(sparse_meshes[i], f'{TEMPLATE_OUTPUT_FOLDER}{MESH_FILENAMES[i]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Template itk::ImageToPointSetFilter<itk::Image<signedshort,2>,itk::PointSet<signedshort,2>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterISS2PSSS2'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<signedshort,2>,itk::PointSet<signedshort,2>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<signedshort,3>,itk::PointSet<signedshort,3>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterISS3PSSS3'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<signedshort,3>,itk::PointSet<signedshort,3>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<signedshort,4>,itk::PointSet<signedshort,4>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterISS4PSSS4'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<signedshort,4>,itk::PointSet<signedshort,4>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<unsignedchar,2>,itk::PointSet<unsignedchar,2>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterIUC2PSUC2'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<unsignedchar,2>,itk::PointSet<unsignedchar,2>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<unsignedchar,3>,itk::PointSet<unsignedchar,3>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterIUC3PSUC3'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<unsignedchar,3>,itk::PointSet<unsignedchar,3>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<unsignedchar,4>,itk::PointSet<unsignedchar,4>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterIUC4PSUC4'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<unsignedchar,4>,itk::PointSet<unsignedchar,4>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<unsignedshort,2>,itk::PointSet<unsignedshort,2>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterIUS2PSUS2'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<unsignedshort,2>,itk::PointSet<unsignedshort,2>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<unsignedshort,3>,itk::PointSet<unsignedshort,3>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterIUS3PSUS3'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<unsignedshort,3>,itk::PointSet<unsignedshort,3>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<unsignedshort,4>,itk::PointSet<unsignedshort,4>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterIUS4PSUS4'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<unsignedshort,4>,itk::PointSet<unsignedshort,4>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<float,2>,itk::PointSet<float,2>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterIF2PSF2'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<float,2>,itk::PointSet<float,2>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<float,3>,itk::PointSet<float,3>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterIF3PSF3'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<float,3>,itk::PointSet<float,3>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<float,4>,itk::PointSet<float,4>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterIF4PSF4'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<float,4>,itk::PointSet<float,4>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<double,2>,itk::PointSet<double,2>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterID2PSD2'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<double,2>,itk::PointSet<double,2>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<double,3>,itk::PointSet<double,3>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterID3PSD3'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<double,3>,itk::PointSet<double,3>>'\n",
      "Template itk::ImageToPointSetFilter<itk::Image<double,4>,itk::PointSet<double,4>>\n",
      " already defined as <class 'itk.itkImageToPointSetFilterPython.itkImageToPointSetFilterID4PSD4'>\n",
      " is redefined as {cl}\n",
      "Warning: template already defined 'itk::ImageToPointSetFilter<itk::Image<double,4>,itk::PointSet<double,4>>'\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "48bf8486ae32422aa2080915a6b44da6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'numberOfComponents': 3, 'd…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualize with itkwidgets\n",
    "view(geometries=[mesh for mesh in sparse_meshes])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Register Sample for Correspondence\n",
    "Our goal is to find correspondence points between the template and each sample mesh. In order to get correspondence, two steps are employed:\n",
    "- First, a copy of the template mesh is registered to the target sample;\n",
    "- Second, each mesh point is repositioned at its nearest neighbor on the target sample.\n",
    "\n",
    "The result of this process is a full collection of meshes having the same number of points and correspondence between each point such that it represents the same approximate feature on each femur. The cells below show an example of the template being registered to a single sample mesh. This process is repeated for every sample in the full pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration: 0 Metric: 0.0\n",
      "Iteration: 0 Metric: 0.0\n",
      "Iteration: 0 Metric: 0.0\n",
      "Iteration: 0 Metric: 0.0\n",
      "Iteration: 0 Metric: 0.0\n",
      "Iteration: 0 Metric: 0.7779061854942507\n",
      "Iteration: 1 Metric: 0.20232000292614374\n",
      "Iteration: 2 Metric: 0.13674457929331374\n",
      "Iteration: 3 Metric: 0.12514542842076454\n",
      "Iteration: 4 Metric: 0.122550291731051\n",
      "Iteration: 5 Metric: 0.12061420415595286\n",
      "Iteration: 6 Metric: 0.11841144247855505\n",
      "Iteration: 7 Metric: 0.11589176266252948\n",
      "Iteration: 8 Metric: 0.11325031330727194\n",
      "Iteration: 9 Metric: 0.11063534863433273\n",
      "Iteration: 10 Metric: 0.10807184441782834\n",
      "Iteration: 11 Metric: 0.10561519395469256\n",
      "Iteration: 12 Metric: 0.10323751692873938\n",
      "Iteration: 13 Metric: 0.10090638903838738\n",
      "Iteration: 14 Metric: 0.09868776219709642\n",
      "Iteration: 15 Metric: 0.09658259228964905\n",
      "Iteration: 16 Metric: 0.09458964512549776\n",
      "Iteration: 17 Metric: 0.09270177175503178\n",
      "Iteration: 18 Metric: 0.09092719465045711\n",
      "Iteration: 19 Metric: 0.08921088449330158\n",
      "Iteration: 20 Metric: 0.08751331057883362\n",
      "Iteration: 21 Metric: 0.08595625062114805\n",
      "Iteration: 22 Metric: 0.08448978369964545\n",
      "Iteration: 23 Metric: 0.08306276567722867\n",
      "Iteration: 24 Metric: 0.08171543662323255\n",
      "Iteration: 25 Metric: 0.08043608000343186\n",
      "Iteration: 26 Metric: 0.07922975383845902\n",
      "Iteration: 27 Metric: 0.07806911113091398\n",
      "Iteration: 28 Metric: 0.0769828791766969\n",
      "Iteration: 29 Metric: 0.07592302414933781\n",
      "Iteration: 30 Metric: 0.07494589189324528\n",
      "Iteration: 31 Metric: 0.07396736060192158\n",
      "Iteration: 32 Metric: 0.07305007692282282\n",
      "Iteration: 33 Metric: 0.07217091810544289\n",
      "Iteration: 34 Metric: 0.07137870580845175\n",
      "Iteration: 35 Metric: 0.07062566206862178\n",
      "Iteration: 36 Metric: 0.06990448486572416\n",
      "Iteration: 37 Metric: 0.06922109639851765\n",
      "Iteration: 38 Metric: 0.06853380709305265\n",
      "Iteration: 39 Metric: 0.06788873436392066\n",
      "Iteration: 40 Metric: 0.0672722282894907\n",
      "Iteration: 41 Metric: 0.06669437953985918\n",
      "Iteration: 42 Metric: 0.06615220704647713\n",
      "Iteration: 43 Metric: 0.06561881914708118\n",
      "Iteration: 44 Metric: 0.0651103655497439\n",
      "Iteration: 45 Metric: 0.06463178159354313\n",
      "Iteration: 46 Metric: 0.06419378162288007\n",
      "Iteration: 47 Metric: 0.06378539931051558\n",
      "Iteration: 48 Metric: 0.06339390676537018\n",
      "Iteration: 49 Metric: 0.06299462478066936\n",
      "Iteration: 50 Metric: 0.06259765910714457\n",
      "Iteration: 51 Metric: 0.062223501410776164\n",
      "Iteration: 52 Metric: 0.06186476948546618\n",
      "Iteration: 53 Metric: 0.06151061722063784\n",
      "Iteration: 54 Metric: 0.06118354489978483\n",
      "Iteration: 55 Metric: 0.06087213131634895\n",
      "Iteration: 56 Metric: 0.06057013719210121\n",
      "Iteration: 57 Metric: 0.06027498100856217\n",
      "Iteration: 58 Metric: 0.05999541178989893\n",
      "Iteration: 59 Metric: 0.05972059722194064\n",
      "Iteration: 60 Metric: 0.05944584050731294\n",
      "Iteration: 61 Metric: 0.05917541864616805\n",
      "Iteration: 62 Metric: 0.058908971747794144\n",
      "Iteration: 63 Metric: 0.058655260566518784\n",
      "Iteration: 64 Metric: 0.058415430581464645\n",
      "Iteration: 65 Metric: 0.058178836264712354\n",
      "Iteration: 66 Metric: 0.057937409843928145\n",
      "Iteration: 67 Metric: 0.057697151558410986\n",
      "Iteration: 68 Metric: 0.05747225797292385\n",
      "Iteration: 69 Metric: 0.05727423976713176\n",
      "Iteration: 70 Metric: 0.05709037866354649\n",
      "Iteration: 71 Metric: 0.05691542628966276\n",
      "Iteration: 72 Metric: 0.056746657780375326\n",
      "Iteration: 73 Metric: 0.05657913987842849\n",
      "Iteration: 74 Metric: 0.05641514413854398\n",
      "Iteration: 75 Metric: 0.056251362421798484\n",
      "Iteration: 76 Metric: 0.056090432428837184\n",
      "Iteration: 77 Metric: 0.05593873693335362\n",
      "Iteration: 78 Metric: 0.05579209223722084\n",
      "Iteration: 79 Metric: 0.05565609855107295\n",
      "Iteration: 80 Metric: 0.0555228740624395\n",
      "Iteration: 81 Metric: 0.05538491139523176\n",
      "Iteration: 82 Metric: 0.055250594503693316\n",
      "Iteration: 83 Metric: 0.055121970213393746\n",
      "Iteration: 84 Metric: 0.05500291647599293\n",
      "Iteration: 85 Metric: 0.054891744576532966\n",
      "Iteration: 86 Metric: 0.05478383310930921\n",
      "Iteration: 87 Metric: 0.054681873426192164\n",
      "Iteration: 88 Metric: 0.0545842727694528\n",
      "Iteration: 89 Metric: 0.054489492014450656\n",
      "Iteration: 90 Metric: 0.05438684625230607\n",
      "Iteration: 91 Metric: 0.05429039237273384\n",
      "Iteration: 92 Metric: 0.05420227931773402\n",
      "Iteration: 93 Metric: 0.0541178113578612\n",
      "Iteration: 94 Metric: 0.05402946310069399\n",
      "Iteration: 95 Metric: 0.0539281294698308\n",
      "Iteration: 96 Metric: 0.05383914528785153\n",
      "Iteration: 97 Metric: 0.053758757355710056\n",
      "Iteration: 98 Metric: 0.053674653015075825\n",
      "Iteration: 99 Metric: 0.053584798945789866\n",
      "Iteration: 100 Metric: 0.05349382055175694\n",
      "Iteration: 101 Metric: 0.05341009320456643\n",
      "Iteration: 102 Metric: 0.05332881383329767\n",
      "Iteration: 103 Metric: 0.053256532548124506\n",
      "Iteration: 104 Metric: 0.05319328769701666\n",
      "Iteration: 105 Metric: 0.053131653243475706\n",
      "Iteration: 106 Metric: 0.05307704756507718\n",
      "Iteration: 107 Metric: 0.05302060581934514\n",
      "Iteration: 108 Metric: 0.05295340984159919\n",
      "Iteration: 109 Metric: 0.05289001870842689\n",
      "Iteration: 110 Metric: 0.05282385749602994\n",
      "Iteration: 111 Metric: 0.052765058546700064\n",
      "Iteration: 112 Metric: 0.052697587776129824\n",
      "Iteration: 113 Metric: 0.05263494878354366\n",
      "Iteration: 114 Metric: 0.05256945628213697\n",
      "Iteration: 115 Metric: 0.05250504946239966\n",
      "Iteration: 116 Metric: 0.052450526747697526\n",
      "Iteration: 117 Metric: 0.05239895469950103\n",
      "Iteration: 118 Metric: 0.05235319450752982\n",
      "Iteration: 119 Metric: 0.052312412068669024\n",
      "Iteration: 120 Metric: 0.052268564087777705\n",
      "Iteration: 121 Metric: 0.0522214106949559\n",
      "Iteration: 122 Metric: 0.052173354704088175\n",
      "Iteration: 123 Metric: 0.05212340589885686\n",
      "Iteration: 124 Metric: 0.052078422799893984\n",
      "Iteration: 125 Metric: 0.052031849859089076\n",
      "Iteration: 126 Metric: 0.051988323846011\n",
      "Iteration: 127 Metric: 0.051949387451224995\n",
      "Iteration: 128 Metric: 0.05191534420597778\n",
      "Iteration: 129 Metric: 0.05188183736711797\n",
      "Iteration: 130 Metric: 0.051841031615754214\n",
      "Iteration: 131 Metric: 0.05180304581125652\n",
      "Iteration: 132 Metric: 0.05177053415449397\n",
      "Iteration: 133 Metric: 0.05174146628941411\n",
      "Iteration: 134 Metric: 0.05170791558239239\n",
      "Iteration: 135 Metric: 0.05168041336473835\n",
      "Iteration: 136 Metric: 0.051653407663689116\n",
      "Iteration: 137 Metric: 0.05162798857126075\n",
      "Iteration: 138 Metric: 0.051605632399284146\n",
      "Iteration: 139 Metric: 0.05158414454641372\n",
      "Iteration: 140 Metric: 0.05156124141667879\n",
      "Iteration: 141 Metric: 0.051535505328572255\n",
      "Iteration: 142 Metric: 0.05151744439212391\n",
      "Iteration: 143 Metric: 0.0514928965718451\n",
      "Iteration: 144 Metric: 0.05146978819894209\n",
      "Iteration: 145 Metric: 0.051452394271285125\n",
      "Iteration: 146 Metric: 0.05143333849263162\n",
      "Iteration: 147 Metric: 0.05141683983858978\n",
      "Iteration: 148 Metric: 0.05140006000618391\n",
      "Iteration: 149 Metric: 0.051379632528070346\n",
      "Iteration: 150 Metric: 0.051364163225486044\n",
      "Iteration: 151 Metric: 0.051344148804672154\n",
      "Iteration: 152 Metric: 0.0513209198248402\n",
      "Iteration: 153 Metric: 0.051300204451203146\n",
      "Iteration: 154 Metric: 0.051279847038341424\n",
      "Iteration: 155 Metric: 0.05125897474166499\n",
      "Iteration: 156 Metric: 0.05124039364195694\n",
      "Iteration: 157 Metric: 0.051225711659785685\n",
      "Iteration: 158 Metric: 0.051213232380085696\n",
      "Iteration: 159 Metric: 0.05119908540041881\n",
      "Iteration: 160 Metric: 0.05118816073223949\n",
      "Iteration: 161 Metric: 0.05117790965210165\n",
      "Iteration: 162 Metric: 0.051164989486057126\n",
      "Iteration: 163 Metric: 0.05115162286003827\n",
      "Iteration: 164 Metric: 0.051141654999141875\n",
      "Iteration: 165 Metric: 0.05112915649279891\n",
      "Iteration: 166 Metric: 0.05112134513222777\n",
      "Iteration: 167 Metric: 0.051113475657624995\n",
      "Iteration: 168 Metric: 0.05110545616598484\n",
      "Iteration: 169 Metric: 0.05109767407846486\n",
      "Iteration: 170 Metric: 0.051086077346439454\n",
      "Iteration: 171 Metric: 0.05107590830530297\n",
      "Iteration: 172 Metric: 0.05106466261954617\n",
      "Iteration: 173 Metric: 0.05105306999550187\n",
      "Iteration: 174 Metric: 0.05104524411816563\n",
      "Iteration: 175 Metric: 0.051037250720081036\n",
      "Iteration: 176 Metric: 0.051027554599183546\n",
      "Iteration: 177 Metric: 0.051019245272718464\n",
      "Iteration: 178 Metric: 0.0510103015813352\n",
      "Iteration: 179 Metric: 0.05099999631723761\n",
      "Iteration: 180 Metric: 0.05099747702422614\n",
      "Iteration: 181 Metric: 0.05099025910556314\n",
      "Iteration: 182 Metric: 0.05098246677270456\n",
      "Iteration: 183 Metric: 0.05097646728829948\n",
      "Iteration: 184 Metric: 0.0509680084125315\n",
      "Iteration: 185 Metric: 0.050957756849324895\n",
      "Iteration: 186 Metric: 0.05094888751493839\n",
      "Iteration: 187 Metric: 0.05093697704947861\n",
      "Iteration: 188 Metric: 0.05092846670001227\n",
      "Iteration: 189 Metric: 0.050923196874734906\n",
      "Iteration: 190 Metric: 0.0509168646827519\n",
      "Iteration: 191 Metric: 0.050905025166283834\n",
      "Iteration: 192 Metric: 0.05088763831250859\n",
      "Iteration: 193 Metric: 0.05088060842698278\n",
      "Iteration: 194 Metric: 0.05087075561956609\n",
      "Iteration: 195 Metric: 0.050859520430219664\n",
      "Iteration: 196 Metric: 0.05085411683245611\n",
      "Iteration: 197 Metric: 0.05084410138767771\n",
      "Iteration: 198 Metric: 0.05083532634158362\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration: 199 Metric: 0.050829382864870494\n",
      "Iteration: 200 Metric: 0.05082170914386962\n",
      "Iteration: 201 Metric: 0.05081423866472783\n",
      "Iteration: 202 Metric: 0.05080568023614978\n",
      "Iteration: 203 Metric: 0.05079902262172653\n",
      "Iteration: 204 Metric: 0.05079059523282873\n",
      "Iteration: 205 Metric: 0.05078300527810398\n",
      "Iteration: 206 Metric: 0.05077665585338589\n",
      "Iteration: 207 Metric: 0.05076873757269327\n",
      "Iteration: 208 Metric: 0.05076176368430392\n",
      "Iteration: 209 Metric: 0.05076191191460448\n",
      "Iteration: 210 Metric: 0.050759255781437584\n",
      "Iteration: 211 Metric: 0.05075297240683363\n",
      "Iteration: 212 Metric: 0.05074280619842757\n",
      "Iteration: 213 Metric: 0.0507353605396458\n",
      "Iteration: 214 Metric: 0.0507293774424416\n",
      "Iteration: 215 Metric: 0.05072597680923839\n",
      "Iteration: 216 Metric: 0.050725369103616885\n",
      "Iteration: 217 Metric: 0.0507235730125808\n",
      "Iteration: 218 Metric: 0.05071929349275945\n",
      "Iteration: 219 Metric: 0.05071276157849241\n",
      "Iteration: 220 Metric: 0.050704833611047065\n",
      "Iteration: 221 Metric: 0.05070007184266171\n",
      "Iteration: 222 Metric: 0.05069326759814266\n",
      "Iteration: 223 Metric: 0.050687986995510546\n",
      "Iteration: 224 Metric: 0.05068121272230936\n",
      "Iteration: 225 Metric: 0.05067712628603081\n",
      "Iteration: 226 Metric: 0.050671607763591114\n",
      "Iteration: 227 Metric: 0.05066573637544863\n",
      "Iteration: 228 Metric: 0.050661630055786895\n",
      "Iteration: 229 Metric: 0.05065762856878695\n",
      "Iteration: 230 Metric: 0.050653153150137\n",
      "Iteration: 231 Metric: 0.05064987299222196\n",
      "Iteration: 232 Metric: 0.05064315389881709\n",
      "Iteration: 233 Metric: 0.05063952433441711\n",
      "Iteration: 234 Metric: 0.05063466540346055\n",
      "Iteration: 235 Metric: 0.05062943284705007\n",
      "Iteration: 236 Metric: 0.05062446513810828\n",
      "Iteration: 237 Metric: 0.0506187349444279\n",
      "Iteration: 238 Metric: 0.05061558707153199\n",
      "Iteration: 239 Metric: 0.05061220664007433\n",
      "Iteration: 240 Metric: 0.050608453663761495\n",
      "Iteration: 241 Metric: 0.05060405388796036\n",
      "Iteration: 242 Metric: 0.05059907738048261\n",
      "Iteration: 243 Metric: 0.05059543121583007\n",
      "Iteration: 244 Metric: 0.05059200253936908\n",
      "Iteration: 245 Metric: 0.05059161578220578\n",
      "Iteration: 246 Metric: 0.05058699150323093\n",
      "Iteration: 247 Metric: 0.05058291307147369\n",
      "Iteration: 248 Metric: 0.050577805145974174\n",
      "Iteration: 249 Metric: 0.05057217165358974\n",
      "Iteration: 250 Metric: 0.050571062925596905\n",
      "Iteration: 251 Metric: 0.05056496573044829\n",
      "Iteration: 252 Metric: 0.050562110970127035\n",
      "Iteration: 253 Metric: 0.05055728536395525\n",
      "Iteration: 254 Metric: 0.05055235366588629\n",
      "Iteration: 255 Metric: 0.05054737396130461\n",
      "Iteration: 256 Metric: 0.05054440043124623\n",
      "Iteration: 257 Metric: 0.05054188303247858\n",
      "Iteration: 258 Metric: 0.05053881837291953\n",
      "Iteration: 259 Metric: 0.05053234625570165\n",
      "Iteration: 260 Metric: 0.05052936539212418\n",
      "Iteration: 261 Metric: 0.05052522606071205\n",
      "Iteration: 262 Metric: 0.050521127631690255\n",
      "Iteration: 263 Metric: 0.05052246227319725\n",
      "Iteration: 264 Metric: 0.050517423605369864\n",
      "Iteration: 265 Metric: 0.05051349205007255\n",
      "Iteration: 266 Metric: 0.05051099265176425\n",
      "Iteration: 267 Metric: 0.050508973692926215\n",
      "Iteration: 268 Metric: 0.05050551462038012\n",
      "Iteration: 269 Metric: 0.050502524037637345\n",
      "Iteration: 270 Metric: 0.050499477928054426\n",
      "Iteration: 271 Metric: 0.05049629849772875\n",
      "Iteration: 272 Metric: 0.05049760118512804\n",
      "Iteration: 273 Metric: 0.050494756776694634\n",
      "Iteration: 274 Metric: 0.05049055434412369\n",
      "Iteration: 275 Metric: 0.05048739284126354\n",
      "Iteration: 276 Metric: 0.050482198188618564\n",
      "Iteration: 277 Metric: 0.05047929580808147\n",
      "Iteration: 278 Metric: 0.0504751878945688\n",
      "Iteration: 279 Metric: 0.05047077602731021\n",
      "Iteration: 280 Metric: 0.050465679778863945\n",
      "Iteration: 281 Metric: 0.050462028373489266\n",
      "Iteration: 282 Metric: 0.05045805846907594\n",
      "Iteration: 283 Metric: 0.050454642707376894\n",
      "Iteration: 284 Metric: 0.050450230167405796\n",
      "Iteration: 285 Metric: 0.050446779213573344\n",
      "Iteration: 286 Metric: 0.050443380421668425\n",
      "Iteration: 287 Metric: 0.05044124380681699\n",
      "Iteration: 288 Metric: 0.0504394860849241\n",
      "Iteration: 289 Metric: 0.05043804286975267\n",
      "Iteration: 290 Metric: 0.050437299696372666\n",
      "Iteration: 291 Metric: 0.050432846525755165\n",
      "Iteration: 292 Metric: 0.05043123811402428\n",
      "Iteration: 293 Metric: 0.05042600614502405\n",
      "Iteration: 294 Metric: 0.05042245290959345\n",
      "Iteration: 295 Metric: 0.050417389038791026\n",
      "Iteration: 296 Metric: 0.05041356076936511\n",
      "Iteration: 297 Metric: 0.05040913396334622\n",
      "Iteration: 298 Metric: 0.05040420966281016\n",
      "Iteration: 299 Metric: 0.05040278224274329\n",
      "Iteration: 300 Metric: 0.050401768630244766\n",
      "Iteration: 301 Metric: 0.05039872658215754\n",
      "Iteration: 302 Metric: 0.05039429292407768\n",
      "Iteration: 303 Metric: 0.05039157619714431\n",
      "Iteration: 304 Metric: 0.05038822246984111\n",
      "Iteration: 305 Metric: 0.050383375177383075\n",
      "Iteration: 306 Metric: 0.0503805094230236\n",
      "Iteration: 307 Metric: 0.05037832060080366\n",
      "Iteration: 308 Metric: 0.050374129488421786\n",
      "Iteration: 309 Metric: 0.05037167465089316\n",
      "Iteration: 310 Metric: 0.050368475933972694\n",
      "Iteration: 311 Metric: 0.05036502987680338\n",
      "Iteration: 312 Metric: 0.05036328948124166\n",
      "Iteration: 313 Metric: 0.05036304455490719\n",
      "Iteration: 314 Metric: 0.050360530034586085\n",
      "Iteration: 315 Metric: 0.05035837642959368\n",
      "Iteration: 316 Metric: 0.05035485754947832\n",
      "Iteration: 317 Metric: 0.05035390608912031\n",
      "Iteration: 318 Metric: 0.05035212109171042\n",
      "Iteration: 319 Metric: 0.05035103871150679\n",
      "Iteration: 320 Metric: 0.05035011328891645\n",
      "Iteration: 321 Metric: 0.05034859021203525\n",
      "Iteration: 322 Metric: 0.050347311419667286\n",
      "Iteration: 323 Metric: 0.05034482451091769\n",
      "Iteration: 324 Metric: 0.050344502999895985\n",
      "Iteration: 325 Metric: 0.05034167345627448\n",
      "Iteration: 326 Metric: 0.05033954454776951\n",
      "Iteration: 327 Metric: 0.050338444807378854\n",
      "Iteration: 328 Metric: 0.05033333892514485\n",
      "Iteration: 329 Metric: 0.05033061695326114\n",
      "Iteration: 330 Metric: 0.050328536986903295\n",
      "Iteration: 331 Metric: 0.050325284955912646\n",
      "Iteration: 332 Metric: 0.05032113795012483\n",
      "Iteration: 333 Metric: 0.050318440363230674\n",
      "Iteration: 334 Metric: 0.050315692835759554\n",
      "Iteration: 335 Metric: 0.050310221774776884\n",
      "Iteration: 336 Metric: 0.050307172757623336\n",
      "Iteration: 337 Metric: 0.050304335918295265\n",
      "Iteration: 338 Metric: 0.050302998852957415\n",
      "Iteration: 339 Metric: 0.050299817591392314\n",
      "Iteration: 340 Metric: 0.05029759969848791\n",
      "Iteration: 341 Metric: 0.05029510993358949\n",
      "Iteration: 342 Metric: 0.05029122136606211\n",
      "Iteration: 343 Metric: 0.05028921697481013\n",
      "Iteration: 344 Metric: 0.05028675933115594\n",
      "Iteration: 345 Metric: 0.050283812442764295\n",
      "Iteration: 346 Metric: 0.050281927743758954\n",
      "Iteration: 347 Metric: 0.05027902561494371\n",
      "Iteration: 348 Metric: 0.05027841284200891\n",
      "Iteration: 349 Metric: 0.050273692579973356\n",
      "Iteration: 350 Metric: 0.05027127734164066\n",
      "Iteration: 351 Metric: 0.05026854623742124\n",
      "Iteration: 352 Metric: 0.05026598477111682\n",
      "Iteration: 353 Metric: 0.05026377775027612\n",
      "Iteration: 354 Metric: 0.05026080899514568\n",
      "Iteration: 355 Metric: 0.050260336182794914\n",
      "Iteration: 356 Metric: 0.05025836559014859\n",
      "Iteration: 357 Metric: 0.050256441618175454\n",
      "Iteration: 358 Metric: 0.05025404058669818\n",
      "Iteration: 359 Metric: 0.050251581225139806\n",
      "Iteration: 360 Metric: 0.050250054709434146\n",
      "Iteration: 361 Metric: 0.05024844423755407\n",
      "Iteration: 362 Metric: 0.05024696256938067\n",
      "Iteration: 363 Metric: 0.050243679326953006\n",
      "Iteration: 364 Metric: 0.05024205871054532\n",
      "Iteration: 365 Metric: 0.05024065722054119\n",
      "Iteration: 366 Metric: 0.05023820875244211\n",
      "Iteration: 367 Metric: 0.05023613769674013\n",
      "Iteration: 368 Metric: 0.05023555529977827\n",
      "Iteration: 369 Metric: 0.050234782085916294\n",
      "Iteration: 370 Metric: 0.05023372316288263\n",
      "Iteration: 371 Metric: 0.05023041006198993\n",
      "Iteration: 372 Metric: 0.0502283742613524\n",
      "Iteration: 373 Metric: 0.05022586075701607\n",
      "Iteration: 374 Metric: 0.050223304015505024\n",
      "Iteration: 375 Metric: 0.05022193669254255\n",
      "Iteration: 376 Metric: 0.05021741490878071\n",
      "Iteration: 377 Metric: 0.05021548423049242\n",
      "Iteration: 378 Metric: 0.05021318832704701\n",
      "Iteration: 379 Metric: 0.05021125341836687\n",
      "Iteration: 380 Metric: 0.050210314308740976\n",
      "Iteration: 381 Metric: 0.050207382476498016\n",
      "Iteration: 382 Metric: 0.050205047603897116\n",
      "Iteration: 383 Metric: 0.05020377646884177\n",
      "Iteration: 384 Metric: 0.050202725803133\n",
      "Iteration: 385 Metric: 0.05020159735330152\n",
      "Iteration: 386 Metric: 0.05019977700818119\n",
      "Iteration: 387 Metric: 0.05019659475499572\n",
      "Iteration: 388 Metric: 0.05019382863798433\n",
      "Iteration: 389 Metric: 0.0501914220347015\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration: 390 Metric: 0.05018866379176612\n",
      "Iteration: 391 Metric: 0.0501866580990529\n",
      "Iteration: 392 Metric: 0.05018544461556436\n",
      "Iteration: 393 Metric: 0.050183944161473114\n",
      "Iteration: 394 Metric: 0.05018252535440755\n",
      "Iteration: 395 Metric: 0.050181952039632355\n",
      "Iteration: 396 Metric: 0.05018052191351881\n",
      "Iteration: 397 Metric: 0.05017906303736656\n",
      "Iteration: 398 Metric: 0.050178842342465696\n",
      "Iteration: 399 Metric: 0.05017926878989274\n",
      "Iteration: 400 Metric: 0.0501761496165269\n",
      "Iteration: 401 Metric: 0.05017516704134721\n",
      "Iteration: 402 Metric: 0.05017371936068065\n",
      "Iteration: 403 Metric: 0.05017193665910303\n",
      "Iteration: 404 Metric: 0.050170784743674424\n",
      "Iteration: 405 Metric: 0.050170514887143504\n",
      "Iteration: 406 Metric: 0.050171765973031336\n",
      "Iteration: 407 Metric: 0.050170869118313775\n",
      "Iteration: 408 Metric: 0.05016937934734629\n",
      "Iteration: 409 Metric: 0.0501673844188967\n",
      "Iteration: 410 Metric: 0.050166100196711874\n",
      "Iteration: 411 Metric: 0.05016514666894262\n",
      "Iteration: 412 Metric: 0.05016612976804814\n",
      "Iteration: 413 Metric: 0.05016357882674013\n",
      "Iteration: 414 Metric: 0.050162055361429175\n",
      "Iteration: 415 Metric: 0.05016021897741215\n",
      "Iteration: 416 Metric: 0.050158974869405\n",
      "Iteration: 417 Metric: 0.050158220964876724\n",
      "Iteration: 418 Metric: 0.050157306725792966\n",
      "Iteration: 419 Metric: 0.050155201441845036\n",
      "Iteration: 420 Metric: 0.050160625074469035\n",
      "Iteration: 421 Metric: 0.05015604999918623\n",
      "Iteration: 422 Metric: 0.05015480240021365\n",
      "Iteration: 423 Metric: 0.05015077919491364\n",
      "Iteration: 424 Metric: 0.05014984714187697\n",
      "Iteration: 425 Metric: 0.050147897342842265\n",
      "Iteration: 426 Metric: 0.05014750205717864\n",
      "Iteration: 427 Metric: 0.05014676747573058\n",
      "Iteration: 428 Metric: 0.05014439829130655\n",
      "Iteration: 429 Metric: 0.05014171150148085\n",
      "Iteration: 430 Metric: 0.05013783652494085\n",
      "Iteration: 431 Metric: 0.050136598736406475\n",
      "Iteration: 432 Metric: 0.050136189352769554\n",
      "Iteration: 433 Metric: 0.05013564703574842\n",
      "Iteration: 434 Metric: 0.050137191495986944\n",
      "Iteration: 435 Metric: 0.050135733190428446\n",
      "Iteration: 436 Metric: 0.05013407258124941\n",
      "Iteration: 437 Metric: 0.050131244740288856\n",
      "Iteration: 438 Metric: 0.05013077696487938\n",
      "Iteration: 439 Metric: 0.05012929740093765\n",
      "Iteration: 440 Metric: 0.0501272224464538\n",
      "Iteration: 441 Metric: 0.05012622208074147\n",
      "Iteration: 442 Metric: 0.05012479969805858\n",
      "Iteration: 443 Metric: 0.050124138938231866\n",
      "Iteration: 444 Metric: 0.0501232227145544\n",
      "Iteration: 445 Metric: 0.050122276680540256\n",
      "Iteration: 446 Metric: 0.05012011797450067\n",
      "Iteration: 447 Metric: 0.05011642153928056\n",
      "Iteration: 448 Metric: 0.05011469844334742\n",
      "Iteration: 449 Metric: 0.05011274779865784\n",
      "Iteration: 450 Metric: 0.050110454865674454\n",
      "Iteration: 451 Metric: 0.05010929522282516\n",
      "Iteration: 452 Metric: 0.05010668016619564\n",
      "Iteration: 453 Metric: 0.05010551075639925\n",
      "Iteration: 454 Metric: 0.050103470218839916\n",
      "Iteration: 455 Metric: 0.05010128360402523\n",
      "Iteration: 456 Metric: 0.05009958402536302\n",
      "Iteration: 457 Metric: 0.05009875797171815\n",
      "Iteration: 458 Metric: 0.050099827762251886\n",
      "Iteration: 459 Metric: 0.050098532874809996\n",
      "Iteration: 460 Metric: 0.0500959213737748\n",
      "Iteration: 461 Metric: 0.050094885953764406\n",
      "Iteration: 462 Metric: 0.05009225552261694\n",
      "Iteration: 463 Metric: 0.05009018220725623\n",
      "Iteration: 464 Metric: 0.050089567358585066\n",
      "Iteration: 465 Metric: 0.05008716836496022\n",
      "Iteration: 466 Metric: 0.05008466078994798\n",
      "Iteration: 467 Metric: 0.0500803773630747\n",
      "Iteration: 468 Metric: 0.050075036843470366\n",
      "Iteration: 469 Metric: 0.050071654360699755\n",
      "Iteration: 470 Metric: 0.050070076813275695\n",
      "Iteration: 471 Metric: 0.05006682416396176\n",
      "Iteration: 472 Metric: 0.05006398214713756\n",
      "Iteration: 473 Metric: 0.050063357436816475\n",
      "Iteration: 474 Metric: 0.050062277529736375\n",
      "Iteration: 475 Metric: 0.05006005400160558\n",
      "Iteration: 476 Metric: 0.050057934210079764\n",
      "Iteration: 477 Metric: 0.050056958310114794\n",
      "Iteration: 478 Metric: 0.05005436748371105\n",
      "Iteration: 479 Metric: 0.05005406283415361\n",
      "Iteration: 480 Metric: 0.05005277249189518\n",
      "Iteration: 481 Metric: 0.050051048384666064\n",
      "Iteration: 482 Metric: 0.05005044755432356\n",
      "Iteration: 483 Metric: 0.05004907424439143\n",
      "Iteration: 484 Metric: 0.05004852043744791\n",
      "Iteration: 485 Metric: 0.05004650240393613\n",
      "Iteration: 486 Metric: 0.05004436681719005\n",
      "Iteration: 487 Metric: 0.05004304873360667\n",
      "Iteration: 488 Metric: 0.050041428661499326\n",
      "Iteration: 489 Metric: 0.05004045188778732\n",
      "Iteration: 490 Metric: 0.050039639222020954\n",
      "Iteration: 491 Metric: 0.05003783421892475\n",
      "Iteration: 492 Metric: 0.05003649428342877\n",
      "Iteration: 493 Metric: 0.050036245476085435\n",
      "Iteration: 494 Metric: 0.05003664703200972\n",
      "Iteration: 495 Metric: 0.05003430333762224\n",
      "Iteration: 496 Metric: 0.05003277583055426\n",
      "Iteration: 497 Metric: 0.05003165744994026\n",
      "Iteration: 498 Metric: 0.05003009200447867\n",
      "Iteration: 499 Metric: 0.050030346055022885\n",
      "Iteration: 500 Metric: 0.05003009200447867\n",
      "Number of iterations run: 500\n",
      "Final metric value: 0.05003009200447867\n",
      "Final transform position: [0.8601696919525987, -0.06117567675444691, 0.1455632756241091, 0.024242047636822365, 0.9355363281830277, -0.03948992331825653, 0.024086151345177262, -0.003011293128523798, 0.9202711414886503, -0.9714223510362355, -0.3219046794482, -0.17710610686041697]\n",
      "Optimizer scales: [54.022498598098714, 21.996100536727774, 18.232899837112292, 54.022498598098714, 21.996100536727774, 18.232899837112292, 54.02249859809877, 21.996100536727774, 18.23289983711267, 1.0000000000000018, 1.0000000000000018, 1.0000000000000018]\n",
      "Optimizer learning rate: 1.0\n"
     ]
    }
   ],
   "source": [
    "registered_template = align.register_template_to_sample(sparse_meshes[1],\n",
    "                                                  dense_meshes[2],\n",
    "                                                  learning_rate=1.0,\n",
    "                                                  minimum_convergence_value=-0.1,\n",
    "                                                  max_iterations=500,\n",
    "                                                  verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4cc752eb19bb4f1daa26feb935e83d38",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'numberOfComponents': 3, 'd…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view(geometries=[registered_template, dense_meshes[2]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Resample Template for Correspondence\n",
    "\n",
    "With the template mesh registered to the sample we can find correspondence points on the surface of the sample bone by getting the nearest sample point to every point on the template mesh surface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e46e2901171c4e22a9a8d8ce2babbfe1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'numberOfComponents': 3, 'd…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sample_correspondence = align.resample_template_from_target(registered_template, dense_meshes[2])\n",
    "view(geometries=[sample_correspondence,dense_meshes[2]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define Procrustes Alignment Parameters\n",
    "Once template meshes have been aligned to represent each input mesh with correspondence we can run Procrustes alignment and get out a mean mesh as the new template. Here we demonstrate aligning two sample meshes as an example, but the full population is used in the atlas generation pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "17d0ebec748b49f5a9cd17480fcb60db",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'numberOfComponents': 3, 'd…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "second_registered_template = align.register_template_to_sample(sparse_meshes[1],\n",
    "                                                               dense_meshes[3],\n",
    "                                                               max_iterations=100)\n",
    "second_sample_correspondence = align.resample_template_from_target(second_registered_template, dense_meshes[3])\n",
    "\n",
    "mean_correspondence = align.get_mean_correspondence_mesh([sample_correspondence,second_sample_correspondence],\n",
    "                                                        convergence_threshold=1.0,\n",
    "                                                        verbose=False)\n",
    "view(geometries=[mean_correspondence])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute Hausdorff Distance\n",
    "\n",
    "We can calculate the Hausdorff distance from the current to previous mesh atlas at each iterative refinement to quantify the amount of change between iterations. In this case the meshes are in correspondence so we get the largest Euclidean distance between any pair of correspondence points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hausdorff distance between initial and updated mesh: 2.018958098608886\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2211d8392c8f4787b1aff5a8769439ce",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'numberOfComponents': 3, 'd…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dist = align.get_pairwise_hausdorff_distance(mean_correspondence,sparse_meshes[1])\n",
    "print(f'Hausdorff distance between initial and updated mesh: {dist}')\n",
    "\n",
    "view(geometries=[sparse_meshes[1],mean_correspondence])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run Iterative Refinement\n",
    "\n",
    "Here we run the previously described registration, deformation, and alignment steps on every sample mesh, iterating the template as we go."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fix the number of iterative atlas refinements.\n",
    "# For the mouse femur data set the atlas tends to converge\n",
    "# after 3-5 iterations.\n",
    "NUM_ITERATIONS = 1\n",
    "\n",
    "# Select index of atlas template to refine.\n",
    "# Initial template should be as close to a \"mean\" shape\n",
    "# as possible, with no distinct surface outliers.\n",
    "TEMPLATE_IDX = 1\n",
    "\n",
    "# Prepare directory to write out atlas iterations.\n",
    "# ex. 'Output/mean/0/901-L-femur-label.obj'\n",
    "for i in range(0,NUM_ITERATIONS):\n",
    "    os.makedirs(MEAN_OUTPUT_FOLDER + str(i) + '\\\\', exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time elapsed for iteration 0: 164.92673254013062\n",
      "Writing to Output/mean//0/902-R-femur-label.obj\n"
     ]
    }
   ],
   "source": [
    "template_mesh = sparse_meshes[TEMPLATE_IDX]\n",
    "for iteration in range(0, NUM_ITERATIONS):\n",
    "    starttime = time.time()\n",
    "    updated_template = align.refine_template_from_population(template_mesh=template_mesh,\n",
    "                                                             target_meshes=dense_meshes,\n",
    "                                                             registration_iterations=200)\n",
    "    endtime = time.time()\n",
    "    print(f'Time elapsed for iteration {iteration}: {endtime - starttime}')\n",
    "    \n",
    "    output_path = f'{MEAN_OUTPUT_FOLDER}/{iteration}/{MESH_FILENAMES[TEMPLATE_IDX]}'\n",
    "    print(f'Writing to {output_path}')\n",
    "    itk.meshwrite(updated_template, output_path)\n",
    "    \n",
    "    template_mesh = updated_template"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final template may be further examined via point distance metrics to determine fitness as an atlas. Here we will simply view the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6145223f116848f794ae1da95a839b5b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'numberOfComponents': 3, 'd…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Visualize \n",
    "view(geometries=[template_mesh])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
